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AN ITERaTIVE METHOD FOR FINDING THE 
TRANSIENT RESPONSE OF AN AUTOMATICALLY 


CONTROLLED MISSILB HAVING NON-LINRAR STIPeNESS 


SUMMARY 


An equation of missile motion in the pitching plane was found, 
assuming a non-linearity in aerodynamic stiffness in the form of a 
cubic term. An iterative method for approximating the angle-of- 
attack response to a step displacement of elevator angle was 
formulated. The results of applications of this method were 
compared with computer solutions; convergence of the approximations 


us 
woue demonstrated. 


A simple feedback control system was designed, and the iterative 


method was extended to the case of a control system step input. 


The method was thought to have limited practical value because 


of the labor required in using it. 
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AN ITERATIVE METHOD #OR FINDING THE 
TRANSIENT RESPONSE OF AN AUTOMATICALLY 


CONTROLLED MISSILE HAVING NON-LINEAR STIFFNESS 


INTRODUCTION 


Important non-linearities which are often present in low aspect 
ratio configurations are those of the lift versus incidence and the 
pitching moment versus incidence characteristics. In missile design 
these non-linearities, as well as any others present, are likely to 
be ignored initially, simply because linear system design is well 
developed and relatively easy. Automatic computing techniques may 
then be used to determine the behaviour of the system when the non- 


linearities are included. 


It would be a great help for the designer to be able to find 
the transient response of a contemplated design without going to the 
computer. It was for this reason that the present investigation 
was undertaken. By expressing the pitching moment as the sum of a 
linear term and a cubic term in angle~of-attack, ignoring the 
non-linearity in lift because it is likely to be less important, 
an equation amenable to iterative solution can be formulated. The 
first approximation to the transient response may be expressed as a 
Fourier series, which must be cubed to perform the first iteration. 
The main effort of the development will be to formulate a practical 
method for performing the cubing process, for the first and 
subsequent iterations, and to systematise the steps necessary for 


carrying out the iterations. 


First, an equation relating angle-of-attack to elevator angle 


will be derived, and the response to a step displacement of elevator 





angle will be found. Later, a simple automatic feedback control 
system will be designed for the purpose of formulating a method 


of finding the response to a control system input. 


The author takes this opportunity to express his gratitude for 
the generous assistance given by Mr. P. A. T. Christopher, who 


instigated and supervised this investigation. 





ENGLISH=LETTER NOMENCLATURE 


coefficient of restoring term in Sec. 3 

coefficient of damping term in Sec. 3 

coefficients of terms in the cube of Fourier series 
HM OeC. < 

gain of frequency response at frequency nWpe 

coefficient of cubic term in Sec. 3 

coefficients in the Fourier series for <i, 

phase shift of frequency response at frequency nw 

coefficients in Sec. 2 

phase angles in the Fourier series for Caine 

1/K., times the transfer function of the compensator— 
amplifier 

lift coefficient 

pitching moment coefficient 

partial derivative of My with respect to x, 

coefficient of vane damping term 

decibels 

a/at 

missile body diameter 

actuating error 


coefficients in trigonometric identities 


phase angles in trigonometric identities 





: 
. 
. 
: 


H(s) 


oe 


Ge 


ria 


= 


M'(w) 
Wid (X ) 
Mex 

M x3 


M 


P(s) 


transfer function used in closed-loop iteration 

¥ (-1) 

moment of inertia of missile about y axis 

moment of inertia of vane about its c.g. 

zero frequency gain of the aerodynamic transfer 
function 

zero frequency gain of compensator-amplifier 

gain of input amplifier 

distance between missile c.g. and vane c.g. 

Lart 

mass of missile 

aerodynamic pitching moment on missile about y axis 

aerodynamic pitching moment on vane about its c.g. 

component of moment due to w 

component of moment due to &% 

coefficient of the linear term in the cubic expression 
for pitching moment 

coefficient of the cubic term in the cubic expression 
for pitching moment 

partial derivative of M with respect to n 

harmonic 


1/K, times aerodynamic transfer function 


angular velocity of missile about the y axis 


step function in Sec. 3 





stagnation pressure in appendices I and VI 

magnitude of step function 

non-dimensional step input 

magnitude of non-dimensional step input 

ratio of non-linear to linear steady-state values 

Laplace transform variable 

Wing area 

time 

time measured from start of second half cycle 

period of square wave 

differential operator in equation for missile 
aerodynamics 

component of missile c.g. velocity in x direction 

differential operator in equation for missile 
aerodynamics 

velocity of missile 

vane transfer function 

component of missile c.g. velocity in z direction 

actuator transfer function 

aerodynamic force component in the z direction 

partial derivative of Z with respect to w 

partial derivative of Z with respect to &% 


partial derivative of Z with respect to 1 





ry iw 


uw 


Ke 


GREEK=-LETTER NOMENCLATURE 


missile angle-of-attack 

demanded missile angle-of-attack 

magnitude of input step function 

reference input 

vane angle-of-—attack 

approximations to 

corrections to approximations 

response of linear system to unit step input 

coefficient of the cubic term 

angle between missile axis and vane axis 

damping ratio 

damping ratio of missile 

damping ratio of vane 

damping ratio of elevator actuator 

elevator deflection angle 

magnitude of step elevator inrut 

demanded elevator angle 

angle between missile axis and space reference 

angle between vane axis and space reference 

time constant in the missile aerodynamics transfer 
funeeLon 

fundamental frequency of souare wave 


natural frequency of linear system 





natural frequency of linear missile 
natural frequency of vane 


natural frequency of elevator actuator 


SUBSCRIPTS NOT SHOWN ELSEWHERE 


evaluated function 

function to be derived from an evaluated function 
linear 

non-linear 


steady-state 





SECTION 1: THE EQUATION OF MISSILE MOTION IN PITCH 


The system of axes and sign conventions to be used are 
wijiastrated in Fig. l. Body axes are used. The missile to be 
considered has a fixed wing and an all-moving tail. It may or may 
not be of cruciform configuration, since pure pitching motion will 


be assumed. 


The rigid body equations of motion in the z direction and about 
the y axis (which passes through the center of mass) are, 


respectively, 

ZS m(w - qu) Ge): 
and 

M= Iq (1-2) 


The same equations would, of course, arise if the general three 
dimensional equations of motion were written down and the confining 


assumptions of the present motion imposed. 


Lhe component in the z direction of the total force on the 


missile is taken as 
Z=Zw+Zn (1-3) 


The simplifying assumption has been made here that the gravitational 
force may be neglected because it is small compared with the 
aerodynamic forces available. The contribution to Z due to q is 
omitted on the assumption that it is small. Its inclusion would 
not change the nature of the analysis, but it would be dropped 
eventually anyway because values for it, for the configuration and 
flight conditions considered, are not available. Also, impliciy in 
Eq. 1-3 is the assumption that the malalignment force is zero, 


1.e. Zis zero when w and 2 are zero. 





The aerodynamic moment is taken as 
anit = 
M=M (w) + My ny (1-4) 


where M' (w) is to be read, "The component of moment due to w." 
This is the non-linear term. The same remarks as those following 
Eq. 1-3 apply to the exclusion of components of moment due to q 


and to malalignment. 
Now we take 
u=V (1-5) 


which is valid for small angles of attack, and we assume that 


Vis constant. Then, for small angles of attack, 


xX aor (1-6) 
and 
a = 2K (1-7) 
V 


Substituting Eqs. 1-5 and 1-6 in Eq. 1-l gives 
Z= m(Vx - qV) (1-8) 


Substituting Eqs. 1-6 and 1-7 in Eq. 1-3 gives 


Z= 2, K+ a (1-9) 
Also we replace M' (w) by M'( «&), where M'(«) is to be read, 
"The component of moment due to «." ‘Thus Eq. 1-4 becomes 
= Mit( << = 
M = M"( )+ My (1-10) 


Combining Eqs. 1-8 and 1-9 gives 
Ze A+ ZH = m(Vex - qV) (1-11) 
Combining Eqs. 1-2 and 1-10 gives 


MN( eX) + Mon = Iq (1-12) 
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Eliminating q between Eqs. l-ll and 1l-l2 gives 
mIV x - IZ - Mvi"(&) = IZ, Rn + nVM, Geis) 


Taking 
7 
M"( XS) = My xX + MU of (i= 14) 
equation 1-15 becomes 
oe . 3 ° 
mIVe<~ IZ. % - VM <- nVis X = 12,4 + VN, 9 (1-15) 


Dividing Eq. 1-15 through by mIV gives 





Z M Ne Z M 
ee x, “ bie {er 
Se ie A (1-16) 
mV £ I mV I 
The constants was S. » kK, C and 3 are now defined 
as follows: 
25, 5 (1-17) 
I 
ee a (1-18) 
amV Wm, 
n> 6aee (1-19) 
Tw, 
6 ee (1-20) 
mVK, Wan 
M3 
A= - (ieecay) 
al 


Making these substitutions in Eq. 1-16 gives 
oe 5 u 3 = ey m 
K+ 2 Siwy KX + WY K eB x Se ne ea (1-22) 


In Ref. 1, E. G. Brown-Edwards calculated the aerodynamic 
coefficients for a typical missile configuration. This missile 
1s sketched in Fig. 2. In Appendix I Brown Edwards' results are 
given and are used as the basis of the calculation of the constants 
in kq. l=22. ‘lhe body diameter was taken as one foot and the flight 


Mach number as three. The calculations for a flight altitude of 
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50,000' are shown in the appendix. The results are tabulated in 
Table I. This table also contains the results for sea level 


flight, the calculations for which are not shown. 


SECTION 2: THE TRANSIENT RESPONSE OF THE MISSILE TO A 
STEP DISPLACEMENT OF KLEVATOR ANGLE 


Equation 1-22 relates the angle-of-attack © to the elevator 
angle 4 . We wish to find the variation of % with time upon 
the application of a step displacement in " , initial values all 
being zero. No analytical solution for this non-linear problem 


exists; an iterative method will be sought. 


Equation ]-22 is, in operational form, 


Oe Pk oD + ww, ) of +B? =k) 9, (1 es D)» (2-1) 


This may be written: 


T(D) K + A x? =U(D)y 8), 
where 

7) cet aD 5 Ooi Wan (2-3) 

u(D) = K, Wa (1+. Z D) (2-4) 


Now an approximation to the solution of Eq. 2-2 is the 


solution e<,(t) of 


T(D) X, = U(D) » (2-5) 
This equation is linear, and could be solved analytically, however 


we will use an approximate method. In Ref. 1, where an attempt is 
made to solve the same non-linear problem, Eq. 2-5 is solved 
analytically. When this is done the ensuing steps become quite 


tedious. We hope to find an advantage in taking another course. 
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We find <,(t) in the form of a Fourier series using the 
approximate method given by Wass and Haymen in Ref. 2 and described 
in Appendix II of this thesis. The method requires that the 
frequency response of Eq. 2-5 be determined. The reader who is 
unfamiliar with frequency response technique will find it explained 
in Section 5, in connection with the design of the control system. 


The Wass, Hayman method gives & (t), 
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approximately, in the forn, 
x(t) = Inf [42+ & > 4 sin (mw, t + 3B) (2-6) 
fl ee ali ™ F n 
Te = lie rere sae 


where We is a known constant, |n| is the magnitude of the step 


input, and the A's and B's are known constants - gains and phase 
shifts of the frequency response. If positive damping exists, 


this response is a damped oscillation, with a steady-state value of 
Oe ae (2-7) 


The true steady-state value ( we may be calculated from Eq. 1-22 


by setting all the terms involving derivatives equal to zero: 
2 S aS 
toy, OS ae C=). pa leenan In | (2-8) 


The solution for ( ae is easy, though not analytic. 


Now if we take, as a new approximation to &(t), 


H(t) © (4s X(t) 
Ina, 


(<) Ao + ee | 2-9) 
ee LB 82 | 


then its steady-state value is 


oe) (A) 5 A, = (&) 


O 


S 


Since ay), has the true steady-state value, there is some 
reason for thinking it is an improved approximation for &(t). 
some further justification is presented in Fig. 3. The curves in 
this figure represent the solutions, by digital computer, of an 


equation essentially the same as (1-22). The equation is discussed 





14 


in detail in Section 53. The solutions for the linear case and for 
several non-zero values of the coefficient in the cubic term are 
plotted. In relation to the present example, the linear solution 
corresponds to x(t), and the non-linear solutions to « (t). 

in Fig. 4, each of the non-linear solutions has been scaled up to 
give steady-state values equal to (oe ~- an inversion of the 
process represented by Eq. 2-9. It may be seen that this scaling 
process apparently brings the non-linear solutions into closer 
agreement with the linear one in a general sense, not just the steady- 


state values. ‘hus x ,(t), a scaled-down version of X(t), will 


be taken as an improved approximation to c(t). For brevity let 
x 
R = ( )g i210) 
nla, 
so that 
X(t) = RX (t) (2-11) 


For later convenience let 


x (t)= A042 > aasin (mudt +B.) (2-12) 
V D mM n 


1.e. the response of the linear system to a unit step input. 


aus, 


x(t) In| sere (t) ea 


Hott) = Ri nlX, (t) (2-14) 


Now if £3 is sufficiently small, a closer approximation 


to e(t) is the solution X..(t) of 


TR) X= UD)y - BX? 


Ls _ 5 5 
U(D) n (2 B® |n\ a (221.5 
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Onpuiiet of g(t) is the solution of 
T(D)&X,, = Eee R? ed ae (2-16) 
aul Pa =e I iy 
then by the principle of linear superposition, 
H(t) = X(t) + %,,(t) (2-17) 


To solve (2-16), a method will be developed for expressing 
x? as a Fourier series. 


U 
A A es, 
SE a (pees BON ries Maes ) 
i ; + + Pm FO + BL 


We = ne Ds wae: lis 


The expanded cubic consists of the cube of each of the original 
terms, plus three times the product of the square of each term and 
each of the other terms, plus six times the product of all 
combinations of three unlike terms. Showing a few terms of the 
expansion and introducing a notation for their coefficients we 


have: 


x? 3 2 
— i ud) eee 
is Bo09 A, + Bog) A, A, sin ( ne + B,) “- 
See, 
pecageee gs oan (Wt af Bey 
oa eee i (Wt 4B.) gin (40.44 2.) 4 32. 
aS al F 1 F 3 
+ ne eae sin (iw, t + B. ) sin (jw,,t + B,) sin (kwWt ” B,) (2-18) 


kKach term is identified by a subscript i, j, k. The coefficients 


As 5k take account of the kind of term (e.g. the multiplier 6 for a 


term which is the product of three unlike terms) and the factors 
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24 2/tr , and i/n. The coefficients a are tabulated in Table II. 


Any product of sine terms can be expressed in the forn, 


E, +5 sin (Wt Ee F,) ro 


4 sin (2 wat + F,) 


ec 


eee ae Sane) ey ees 


au 


For example, in the 1,3,5 term we have the identity: 


sin (wt r B, ) sin (Zw ,t a. B,) sin (5w,,t a BS) 


=- 7 sin (Wt -B,-B +B, ) + y sin (3.4 + B,-B_+B 


ahaa 1785 5) 
++ sin (7. -B,+B,+B,) - y sin (9wWt +B, +B,4B, ) 


In this case 


B=-4@ » Esq, E =F @, Eg = ~4 
By = E, = Ey = ES = Be = Ee = Bo — ee eae — ane 0, 
ee ee ane aes 
F_ = a 
3 B B + B. 
Fo = = BL aa B, + Be 
F. = B, +Bo+B 


ge =A 2 2) 


Using this notation a large number of these identities may be 
tabulated concisely, and this has been done, for the combinations 


necessary here, in Table III. 


we adopt the notations: 
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ee rs , iS 
sin (iw,t +B.) sin (j wit + B,) sin (kwit + B,) 


Also, we will write, for later convenience, 


Si | LeU B.) sin (jw, t + B,) 


EF 
Bo(0igj) * =1(0a3) SB (wot + Fa (013)? +... (2-20) 
and 
sin (iwt + B.) 
E (004) & Ey (003) sin (wt i. Fi (coi)? es (2=21) 


The latter expression is a trivial use of such a notation since the 
right hand side reduces to sin (i wt + Be), but its use helps 


provide a systematic tabulation. 


Using the identities represented by Eqs. 2-19, 2-20, and 2-21 


in Kkq. 2-18 gives cay ? as a series of constants and simple sine 
U 


terms. Since sine terms of like frequency can be combined into a 


Single sine term of that frequency, we have: 
Nee, a us = 
i bo + b, sin (wt + c) + b, sin (2 nt + oP) (2-22 ) 
+ woe + De. sin (mw_t + Ce) toes 


B 


where 
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A ? A e A aia. ee 


coe Cole. 


r a ik peGa) A. A. AL. Se 


: 2 
bd, sin (wt + C,) = 80) £1 (001) A, Ay sin (wt + Fy (oo)? 
a E A “A. Sin (wt +F N dare eee 
003 “1(003) “o “3 F 1(003) 
+ a4 sic Ey (45k) A; A, A, sin (w, Xt + Fy (i5k)? oy ee 
b_ sin (n wyt +c et eis 
oe «@ eo ° * ] ww) { * » e@@e@8@ 
+ 8s 5 BA (igk) A; A, A, sin (n Wet + B (agk)? + 
It is possible to shorten the notation conveniently by 
letting 
@“n(ijk) = age Fn(igx) (2-25) 


The values of Cn (ijk) may be determined from the tabulations of 
rake and Ea(isk)* Accordingly, values of c 


are tabulated in Table IV. The coefficients are all listed as 


n(igk) (or simply c) 
positive. ‘The coefficients for which the values of E were negative 
were changed in sign by adding 180° to the corresponding angles F, 
which are also tabulated. It will be seen that like harmonics 


have been grouped together. 


Table IV does not include coefficients for all possible terms. 
The elimination of less important terms was based on the effect of 
each tern on &(t) (the approximation towards which we are working) 
in the numerical example which follows. It had been hoped, at the 
outset of this investigation, that the number of terms needed would 


be small. But it was found necessary to retain a large number of 





ne, 


terms, up to and including 9th harmonic terms. The relative 
importance of the smallest terms retained can best be seen in the 


example of this section. 


The fact that the requirements of a particular numerical 
example must be injected at this point is a deficiency in the method, 
since the choice of retained terms is not necessarily the best choice 
for other examples to which this method will be applied. However, 


to follow a more general course would be impractical. 


To obtain values for oD and Cc for a particular problen, 


the products Cn(isk) A. A. A, and the angles se are 


ies te ijk) 
calculated using the values of gain A and phase shift B from the 
frequency response curves, then like harmonics are combined. The 
combination of several sine terms in the same frequency may be done 
most easily by graphical vector addition. This is described in 


Appendix III. 


The Fourier series for oe may be computed with minimum 
effort by using the form of Table IV. All the necessary results 
of the preceding development are summarised there, and the required 
computations are indicated by its layout. To use the table, 


proceed as follows: 


es Enter the values of gain A and phase shift B of the 
appropriate frequency response. Once this is done, no 
further reference need be made to any information other 


than that contained in Table IV. 
ee Calculate the products A. a Ay 


3. Calculate the products Cn isk) A. A. AL (In the 


table cAAA is written for oe A, AL) 
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a Calculate the angles aie i 


5. Combine like harmonics (e.g. by the method of 
Appendix III), and enter the results, b, and C, , 


in the space provided. 


Step 5 above gives the Fourier series for sate Multiplying 
each of the coefficients of the series by the constant - R? Int? 
gives the Fourier series for - ? In| ? oe ce By Eq. 2-16, 

Xx y(t) may be found by applying the frequency response B/T(aw). 


The next approximation to <(t), if it is desired to go 


further, is X(t), the solution of 


T(a)~, =U(D), -f ee (2-24) 


The determination of a Fourier series noms may be obviated 
as follows. Adding and subtracting pe? on the right hand 
side of (2-24): 


T(D)X, = UD) x ~ px, + 2 (x, - x<.”) (2-25 ) 


let & gx(t) be the solution of 


2 Ss i; 
Then by superposition of the solutions of Eqs. 2-15 and 2-26, 
x ,(t) - X(t) + & (t) (2-27) 


The Fourier series for (<7 oe) need not be determined with 
i) 


as great an accuracy as would be required for that of e To 
solve Eq. 2-26, Xx. and a are each evaluated and cubed, point 


by point. The difference is then plotted, and a numerical harmonic 
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analysis is performed. The frequency response 3/P(ivw) is applied 

to the resulting sourier series, giving & gx: ‘Instead of 

evaluating & a\t) a cyubing Sent pee the series which has been ect 
obtained for we Ber e serie - eee ape a few eee 
is much easier ee. evaluating a Fourler series. But the latter 

method has a certain superiority in that if the series obtained for 

x Ee is not a good representation of the true cube of =; due 

either to a mistake or to a deficiency in the cubing method, then 

the error will be compensated only if the latter method is used. 

A deficiency in the cubing method would probably exist if the linear 


system frequency response was very different from that used in 


formulating the cubing method. 


Improved approximations, ae Ae, etc., would be found in the 


Same way as is oe 


A method is shown in Section 6, in connection with the control 
system response, for simplifying the calculations when the responses 
to more than one input magnitude are required. An analogous 
method can and should be applied here if more than one input 


magnitude is used. 


We now apply this method for the case of missile flight at sea 
level and an input magnitude of .25 radian. This case was chosen 
Sa analog computer solution was available. The constants 
nel Eqs.,c~); and 2-4 are those tabulated for sea level in Table I, 
with one exception: ‘@ is taken as zero. ‘The effect of that 
term on the transient response was found to be negligible, having 
an influence on the frequency response only at very high frequencies, 


and it was omitted from the analog computer setup. It was 


retained, in the development just shown, for generality. 
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The frequency response of kq. 2-5, 


2 
K. Ww 


q (;w)? 4. 243, Ww ie) ee Ww < 
1 i 


= 1 (2-28 ) 
+ --__—_-__ 

i: Ww i. ee yy toa 

Uy 

i 


is plotted in Fig. 5. Gains are plotted in decibels, and phase 
shifts in degrees. Following the Wass, Hayman method for determining 
linear transient response, we take UW), to be one fifth of the 


F 
frequency at peak gain. The peak occurs at 


Ud = ae 


So that 


Oe = 2 Wy 


The gains AL and phase shifts BL of Eq. 2-28, at frequencies 
nw, where n= 0, 1, 2, ..., 13, are tabulated in Columns 2 and 
3 of Table V, with values of n in Column 1. § The values of gain 
are absolute, not decibels. This procedure is necessary to make 
the computations systematic (this differs slightly from the Wass, 
Hayman procedure). 
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The coefficients of the sine terms in the Fourier series for 
a square wave, magnitude .25, are tabulated in Column 4 of Table V. 
Bach coefficient is identified with the harmonic to which it applies 
by the value of n in the same row, e.g. the coefficient of the 
sin 7 Wet term, or seventh harmonic, is .0228. This procedure will 
be followed throughout the thesis, including the use of a second 
column for phase angle when required, and will be referred to as 


"tabulating the Fourier series." 


Having tabulated the Fourier series for the input square wave, 
and the frequency response of the linear system, we find the 
coefficients of the Fourier series for the output, ots by 
multiplying each input coefficient by the corresponding gain. 
the phase angles of the output series are, in this case, the phase 
shifts of the frequency response. The Fourier series for oc) 
is tabulated in Columns 5 and 6. From this tabulation we see, for 
instance, that the fifth harmonic term of the Fourier series for 


X(t) is .0775 sin (5w.t - 88.1"). 


The evaluation of XK, (t), for 39 values of time up to 95 
percent of the half-period of the square wave, is tabulated in 
Table VI. The evaluation was made, as were all of those in this 
report, using the Wass, Hayman method. Seles is plotted in 
Fig. 6 together with the corresponding solution obtained by 
Brown-Edwards on a Short Analog Computer (see Ref. 1). ‘The 
oscillation is seen to be of high amplitude, due to the low damping. 
It is perhaps surprising that the approximation is as good as it is, 
because the Wass, Hayman method assumes that the oscillation will be 
nearly completely damped out by the end of the half cycle of the 


square wave, and it clearly is not. 


Solving Eq. 2-8 for In| = .25 we find 


OE = .195 
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Byetq. 200, 


ae eS Sey, 
By Eq. 2-11, the evaluation of x ,(t) is found from the evaluation 
of x, (t). This is tabulated in Table VI, and plotted, with the 
analog computer solution of the non-linear equation, in Fig. 7. 

x le is calculated by the method of Table IV. This calculation 
is shown in Table VII, and the resulting Fourier series has been 
transferred to Columns 7 and 8 of Table V. Multiplying the 
magnitudes by 


Point oromGoln)? (aes), 


~ 01206 


and adding 180° to the phase angles gives the series for 

- Rn? x, ie or - X,” , and this is tabulated in Columns 9 
and 10. The frequency response ?/T(iw), which in this case only 
differs from U(i wW)/T(G Ww) by a constant multiplier, is tabulated 


in Columns 11 and 12. By Eq. 2-16, “%., is found by applying the 


* 
frequency response (7/T(iwW) to the rate series for - aS , 
and the resulting Fourier series is tabulated in Columns 13 and 14. 
This is evaluated for 27 values of time in Table VI, and, by 

Eq. 2-17, combined with ee to give X(t). ~< (t) is 


plotted in Fig. 7. 


In Fig. 7 we see the relationship between the "reduced" linear 
solution a. the result of the first iteration A(t), and the 
computer solution e<(t). Superficially Ar doesn't look much 
better than & 5» particularly with respect to the amplitudes of 
the peaks and troughs. The times at which the overshoots occur 
have been brought into line however, and it will be shown in Section 3 
that the next iteration gives rapid convergence. The solution to 


the present example was not taken further because the digital computer 
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solution of a similar equation became available to the author, and 
the analog computer solution was thought not to have sufficient 
accuracy to justify its use as a standard with which to compare 


more refined approximations. 






SECTION COMPARISON OF RESULTS WITH SOLUTION OBTAINED WITH 








A DIGITAL COMPUTER 


Mr. P. A. T. Christopher obtained solutions, from a Ferranti 


Mercury digital computer, of the equation: 


Igoe a,Dx + ax + : =o Cay 
where a, = ee. 
a, = 1, 


Q@ is a step function, 


Initial values are zero. 


These solutions were obtained for various valves of b and magnitudes 
of Q. The accuracy of the solutions was such that, for present 


purposes, they may be considered exact. 


The output x(t), of Eq. 3-1, is dependent on four parameters: 
Ay» aos b, and [Ql , the maenitude of the step input. It is shown 
in Avpendix IV that by linear transformations of the variables x 
and t , the output can be made a function of only two parameters: 
the damping ratio § , and the non-dimensional input magnitude IQ] : 
Thus the solution of Eq. 1-22, which is of the same form when 
C@ = 0, could be found from the solution to Eq. 3-1, if ¥ and 
lo for each of the equations were identical. Now the solutions 
to Eq. 3-1 include a range of values of Q.| , but only one value 
of § , i.e. .3, while the value of 5 for the missile at sea level 
is .175. Accepting this deficiency, we will apply the approximation 


method to Eq. 3-1 for the computer value of \c.\ nearest that of 
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the case considered in Section 2. The results will not be 
transformed to the amplitude and time scales of the missile problen, 
because having already seen the approximate missile solution in 
section 2, we are now only interested in the convergence of the 


approximations ~- which does not depend upon scale. 


The equivalent symbols in Eqs. 3-1 and 1-22 are tabulated 
below ( @ = 0): 


Eq. 3-1 Pose 22 
UW) 
a, 2 55 an 
2 
are Pons 
b p 
Z 
Q Ky ce y 


Thus the non-dimensional input in terms relating to Eq. 1-22 is: 


2 

ge SE 5 Meee « Blan 
a 2\5 Ud 

O (w 2 mm 

ie ; 

The value of Ic. corresponding to the case considered in Section 2, 
where the constants are those given in Table I for sea level flight 
Sindee WVn\P = 22255) ce 


.853 i646) 25 


alee aes 
Sol (24.4) 


The computer solutions nearest this are: 


(1) »v = 5 
Voll a> sea 
lol = .4%5 = .283 


(12/2 
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(2) b = 1 
(Oly =. 4 
IQ) = 


van =e OU 
(1)7%e 


We choose the higher value, case ‘Gas and proceed with the application 
of the approximation method to this case. In order to use the 
equations of Section 2 unaltered, we put Eq. 3-1 in the form of 


Eg. 2-2 by letting 


Ci «= x 
fp = be=z1 
T(D) = pe ~ a,D +a = po + .60D + 1 
U(D) = al = 

eee 
" Bie ae 

0 

In| = lel= .4 


The gain of the frequency response 


U(iw 1 
=o. oe (5-2) 
T(iw ) era) ene 


has a peak at W=1, so we take Ud, = .2. Values of gain and 

phase shift of (3-2) are tabulated in Columns 2 and 3 of Table VIII, 
for the required values of NU), The Fourier series for a square 
wave, magnitude .4, is tabulated in Column 4. By Ea. 2-5 we 

find the Fourier series for A(t), and this is tabulated in Columns 

5 and 6. The evaluation of x, (t) for 22 values of time is 

tabulated in Table IX, and is plotted in Fig. 8 with the corresponding 
computer solution. The calculated solution x, appears to be a 


good approximation to the true output of tne linear systen. 
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The steady~state output is found from kq. 2-2 by setting the 


derivative terms equal to zero: 
> 
Qe alii 
For {ny | = 9,4 vunis equation yields 
(X), = .355 
By Eq. 2-10, 


R = = — 888 
a at 
The evaluation of Ce) is found from the evaluation of Soe) 
according to Eq. 2-ll. This is tabulated in Table IX, and is 
plotted in Fig. 9. 


7 


The calculation of ox? is performed in Table X. The 
U 
resulting fourier series is entered in Columns 7 and 8 of Table VIII. 


Multiplying the magnitudes by 


R? Int? = (.888)? (.4)? = .0447 


and adding 180° to the phase angles gives tne series for 
=2° | | , aii or - &<? , and this is tabulated in Columns 9 
and 10. The frequency response #/T(iw) is the same as 
U(iw)/T(iw) in this case because B=U(iw) =1. ‘his has 
already been tabulated in Columns 2 and 35. Thus <5 is found 
according to Eq. 2-16, tabulated in Columns 11 and 12, and 
evaluated in ‘able IX. “x. is obtained using Eq. 2-17, and plotted 
Pry Sig 

The next approximation requires a numerical harmonic analysis. 
In Vable 1X, > xX 2 and x, are cubed, and the difference ( x s - Xx 
is found for each of the 22 values of time. ce ~ Xe | sis 


plotted in Fig. 10. 


3 
3) 





eo 


Now harmonic analysis requires that the function being analysed 
be known throughout one period of the fundamental frequency Wa: 
so far only the first half-cycles have veep Pilar sed ae The 
evaluations of aa and of Jare given, for 22 values of time, in 
Table Al. In this tabulation, the time t, is the time measured 
from the start of the second half—cycle. x<,”, x, and 
( x, - x?) are also tabulated. We see that the values of 
( a Een 3? ) in the second half-cycle are negligible compared 
with those in the first half-cycle, and set them all equal to zero 
in the numerical harmonic analysis. This result is not surprising, 
because the oscillation in the second half-cycle takes place about 
the steady-state value of zero. Therefore the effect of the 
non-linearity is small, and the oscillation of the non-linear system 
is nearly the same as that of the linear one. ‘This can probably be 
taken as a general result, so that ( , _ o..”) need not be 


evaluated for the second half-cycle. 


A 24 ordinate harmonic analysis of ( x,” ~ x7), performed 
using the layout recommended by Wylie in Ref. 3, is shown in 
Table XII. The positions at which the ordinates were taken is 
shown by the vertical lines in Fig. 10. They were found by 
dividing the first half-cycle into 12 equal time intervals. 
The ordinates are labeled Yy9 Yo: etc., in accordance with Wylie's 
notation. The ordinates were all multiplied by 1000 to make the 
harmonic analysis more manageable. The result of the analysis is 
a series in both sine and cosine terms which includes the 12th 
harmonic. Retaining only up to the 9th harmonic, the series was 
put into the form containing only sine terms, and tnis is tabulated 
in Columns 13 and 14 of Table VIII. Kx is obtained by Eq. 2-26, 
and this is tabulated in Columns 15 and 16. It is evaluated in 
Table IX and, in accordance with Eq. 2-27, added to ©“. to give 


3 


Xi This is plotted in Fig. 9 with the other solutions. We 
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see that the convergence is virtually complete in the region of the 
first peak, and probably good enough for most practical purposes in 


the entire range of the computer solution. 


this completes the demonstration of the method develoved in 
section 2. It would be desirable to find out how many iterations 
are necessary for larger values of input magnitude, and for what 
values of input magnitude the first iteration is sufficient, but 


this could not be done by the present author in the time available. 


SECTION 4: EQUATIONS OF COMPONGNTS OF THE MISSILE CONTROL SYSTEM 





A. Incidence Vane 


One of the components of the control system, which will be 
designed in Section 5, 1S an incidence measuring vane at the nose 
of the missile. The design of a conical type vane suitable for 
this use is shown in Appendix V. The vane will be pivoted at its 


Gee.) POSLtion. 


The relationship between the angle of attack of the missile 
and the output of the vane must be determined. Fig. 11 shows the 
angles and sign conventions used in the analysis. the angle Y is 


measured by a transducer in the vane and is the output of the vane. 


It will be seen from the figure that no distinction is made 
between the velocity of the c.g. of the missile and the velocity of 
the vane. This is not strictly correct, of course, because the 
pitching velocity of the missile q creates a component of vane 
velocity perpendicular to the missile axis with magnitude Ly a 
It is easily shown, however, that for the configuration and flight 
conditions under consideration, this component of vane velocity has 
negligible effect on the vane velocity. kor the example being 
considered, under steady-state conditions at an angle of attack of 


40°, the component of vane velocity due to pitching of the missile 
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eS?) 5) ft/sec., and the difference between vare incidence and missile 


Pe : O 
incidence is .l2 . 


The equation of angular motion of the vane about its c.g. is: 


Ms ly Oy (4-1) 


The aerodynamic moment on the vane is taken as 





EN ee en ts (4-2) 
where 
ee ot, 
1 
ae 


and C., sf is the internal damping moment arising from the 
relative ancular motion between the vane and the missile. Ties 


assumed that a device within the vane exists to provide this damping. 


As in the case of the missile (see remarks following Eq. 1-3) 
the aerodynamic damping moment, that due to o<, is omitted on 
the assumption that it is small, but for the real reason that a 


value for it is not available. 


Combining Kgs. 4-1 and 4-2 gives: 


o = S 2 
Cc, X, +c, X 1 (4-3) 
we wish to reduce Eq. 4-3 to a linear equation in the two 

variables < and XX. From the figure, 


o< 


S (4-4) 


R 


and 


Tm eo. = ex (4-5) 
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Substituting (4-4) and (4-5) in (4-3) ives 


¢, (o- ¥) + oe = 1, (8, - x ) (4-6) 
Now = or q represents the pitching accelJeration of the missile, 
and YX that of the vane relative to the missile. It is reasonable 
to assume, and it will be evident later, that the frequency of 
oscillation of the vane is much higher then that of the missile. 
This being the case, and providing that the magnitude of the missile 
oscillation is not excessive, then 6. may be assumed negligible 


in comparison with Nee Accordingly, Eq. 4-6 becomes 


C, (<- ¥)+ oy = -IYX 
_ iy, TORS e C C 
Yr 7p Y- cee a (4-7) 


The constants W, and 5 > are now defined as follows: 
Es 


2 7 VY 


C./ aN Way 


e 
1 


WN 
NO 
i 


iq. 4-/ then becomes 


o@ @ 


¥ + 2 e Sys ree ra ee (4-8) 


The calculation of the constants will be found in Appendix VI. 


They are tabulated in Table I. 
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Be Actuator 


The equation for a typical elevator actuator and the constants 
of that equation were obtained from D. Morton, a student at the 
College of Aeronautics who had investigated the subject. ‘The 
equation is: 


oe 


oes 2 6e Wm Vt em (4-9) 


e 


mm UD 


where Uy is the demanded elevator angle, 
uy is the elevator angle. 


The constants WwW, and 35 are given in Table I. 
2 


SECTION 5: DESIGN OF THE CLOSED-LOOP CONTROL SYSTEM 


A. The Linear System in Transfer Function Form 


The control system arrangement to be used is shown in Fig. l2 


in block diagram form. In this diagram: 


Xp is the demanded missile angle of attack, the control 
system input, 
%g is the reference input, equal to K, times op where 
the constant SS is the gain of the linear amplifier, 
«x is the missile incidence, 
Ny 1s the demanded elevator angle, 
n is the elevator angle, 


¥ is the output of the incidence measuring vane, 


and 


Be er cna (5-1) 


Equation 4-8 relates ™ and Y, Eq. 4-9 relates y and Lp ; 
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and Eq. 1-22 relates o{ and n : Dropping the non-linearity for 
the present so that the design problem is that of a linear system, 


Eq. 1-22 becomes 


Be 5) Wm & + Wns - Kom | n+ 74) (5-2) 


The relation between k and Xp is, as yet, unspecified. 


The Laplace transformation, defined by: 


f(s) = [sit) ro ae 


where f(s) is the epics transform of f(t), will now be applied 
to Eqs. 4-8, 4-9, 5-1, and 5-2, assuming that initial values of 
all the variables and their derivatives are zero. ‘Taking the 
Laplace transforms of these expressions - when initial values are 
zero — simply amounts to substituting s for d/dt and af hor 


a-/at*. These equations then become, after rearranging, 








¥/ rs 
Bie) = vV(4) = Eee PAE ae ict eee (5-3) 
fi i 
am ad 
n(4) Wan 
ee se ee ee (5-4) 
ny (2) az + 2 Dt ae + WwW 
=p 5 % 
E(#) = kK, Se) - ¥ (4) (5-5) 
22 > 
<(e) 5 ple) = 1m, too) (5-6) 
“| (+) 2 ee 


A +2 oo ae Von 
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where V(s), w(s), and P(s) are defined as indicated. whe 


transform of the equation relating E and Ay is given the 





form: 
N p(s) i C(s) i 
BS) K, (5-7) 
where 
ClO)~ Jaf 


In servomechanism terminology, V(s) is the transfer function of 
the vane. Similarly w(s), Ky Pisae and iy Cis) are the 
transfer functions of the actuator, the missile aerodynamics, and 


the compensating network~amplifier. 


Following standard servomechanism technique, the control system 


may be represented by the block diagram in Fig. 13. 


From Eqs. 5-3, 5-4, 5-6, and 5-7 (or directly from Fig. 13): 


(s) = K, K, C(s) W(s) P(s) V(s) (5-3) 


tHtle< 


This is known as the open-loop transfer function, and the system 
comprised of the four elements whose transfer functions appear in it 
and which has an input E(t) and an output Y(t) is the 


"open-loop system". 


Combining Eqs. 5-3, 5-5, and 5-8 to eliminate B(s) and Y(s) 


gives the closed-loop transfer function: 


= K, Ky K, C(s) W(s) P(s) 
Ee a J +K, K, C(s) W(s) P(s) V(s) (5-9) 


C(s) has not been specified. It has been included in order 
to provide the system with satisfactory stability characteristics. 


The stability of the system and the required form for C(s) will 
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be determined using frequency response technique. The choice of 


K, and “ Will first be discussed. 


iB. Steady-State Accuracy 


Suppose that ae is set at some steady value (aes) ee The 
rules of Laplace transforms tell us that (in this type of system) 
o{ takes on a steady value (eX). after the transient motion 
has disappeared, and the ratio oe is given by setting 
s = O in & (s)/%)(s). Accordingly, from Eq. 5-9: 


(x) = K. K 
con S =) a ee ae (5-10) 
S D 1 + K K, 


since 


Ideally the ratio of steady-state values should be one. 
That is the steady-state output should ecual the steady-state 
demand. We see from Hq. 5-10 that this can be accomplished by 
a suitable choice of Ro. regardless of the value of Ky K, : 
But there is another important consideration: if the missile were 
to operate entirely at one flight condition, that is at a particular 
weight, c.g. position, speed, and altitude, then Ky would always 


have the same value. But since missiles rarely operate under such 


conditions, we may assume that Ay will not maintain a sin le 


value throughout the flight, and that (Ge ates will therefore 
vary (unless K, is made to automatically compensate Tess It is 


e, 


evident in Kq. 5-10 that the higher the value of hy Ky can be 
made, the less will be the variation in eu for a given 


K 


ps 


change in K is determined by the aerodynamics of the 


i 1 
missile, but K, may be chosen. 


Aside from considerations of amplifier design, there are two 
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factors which impose an upper limit on K, K One is that as 


1 "2° 


Ky K, is increased, the system becomes less stable, and finally 


unstable. So that the choice of K, must be a compromise between 
stability and steady-state accuracy. The other limit imposed upon 


K, is due to the presence of spurious inputs, or noise, which 


must not be amplified excessively. 


When K. has been chosen, K 


5 may be chosen to give 


y 


x< 
[Jo KK 7 
(Xp) ¢ ee 


for a given flight condition. 


If the steady-state equation for the aerodynamics, including 
the non-linear term, is combined with the steady-state equations 
for the other system components to eliminate all the variables 


except Coe and ey the following equation is obtained: 


/3 : K, K 
Re Fe ap a Ae 
Wy (1+ K, K, ) 1+K, K 


We see that the larger Ky K, is made, the smaller is the 


2 
coefficient of the cubic term. Thus a high open loop gain, seen 
to be desirable in the linear system, reduces the effective 


non-linearity in the non-linear system. 


Ca Ihe Frequency Response of the Linear System 


The frequency response of a linear system at the frequency 
WwW is defined as the complex number obtained by substituting sw for 
Ss in the transfer function relating the output to the input. di 
is usually expressed in polar form, with the magnitude in units of 
decibels. A number is expressed in decibels by multiplying its 


common logarithm by 20. 
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Now in a linear system a steady sinusoidal input of constant 
frequency Ww will give rise to a steady sinusoidal output with 
the same frequency, after the transient motion has disappeared. 
The rules of Laplace transforms tell us that the ratio of the 
output sinusoidal wave amplitude to the input sinusoidal wave 
amplitude, or gain, is equal to the magnitude of the frequency 
response of the system at frequency w; and that the angle by 
which the output sinusoidal wave leads the input sinusoidal wave, 
or phase shift, 1s equal to the argument of the frequency response 
at eeanency w). Furthermore, by virtue of the properties of 
complex numbers, it can be shown that if the transfer function of 
a system is the product of two or more individual transfer functions, 
then the phase shift of the system is equal to the sum of the phase 
shifts of the individual transfer functions; and, if the gains are 
expressed in decibels, then the gain of the system is equal to the 


sum of the gains of the individual transfer functions. 


Ihe frequency response curves of a system, sometimes called 
simply the frequency response, are graphs of the system gain and 


phase shift versus w (or w times some constant ). 


It will be seen that the frequency response of a system may be 
found by substituting iw for D in the differential equation of 
the system, since the transfer function is found by substituting 
s for D and the frequency response by substituting 1W for s. 


This procedure is used elsewhere in this thesis. 


As is shown in textbooks on servomechanisms (see Bibliography), 
the frequency response curves of the open-loop system determine the 
stability of the closed-loop system, and they may be used as the 


basis of trial-and-error design. 


Evaluating the individual freouency response curves and 


combining them in accordance with the rules previously stated, 
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E(iw) V(iw) W(iw) was found and is plotted in Fig. 14. The 
parameter /War is used as the abscissa, and is plotted ona 
logarithmic scale for convenience. The constants are those 
given in Table I for flight at 50,000'. It can be seen from 
Eq. 5-8 that P(iw) V(iw) W(iw) need only be multiplied by 
K K, C(iw) to give the open-loop frequency response. itis 


1 
now necessary to choose a particular compensating network. 


In ordinary servo design the choice of the compensator is 
based on the two main requirements: 1) System stability, and 
2) Independence of the closed-loop system from changes in the 
constants of the components, i.e. high open-loop gain. But in 
the present example an additional requirement must be met: the 
system must have a frequency spread amenable to the Wass, Hayman 
method of transient response calculation (Appendix HE in face, 
it was at this point in the analysis that the writer first became 
aware of the frequency spread requirement, by having chosen a 
compensator that did not satisfy it. To demonstrate: This was 
a lead-lag network with the transfer function: 


(1 + = $)(1 + 


or ieiges, 
C(s) = me | A 

oe de 
Using the logmodang plot, Nichols chart technique (see Bibliography), 
the open-loop gain kK Ky was set at 17 db, giving a closed-loop 
peak gain in the upper range of acceptable values - i.e. a lightly 
damped transient. Be was set at 1.14 db to make the closed-loop 
steady-state gain one. The closed-loop frequency response curves 
are plotted in Fig. 15. The Wass, Hayman method was applied to 
this linear system to find the transient response for a unit step 


input of Xn The fundamental frequency of the input square wave 


was taken as .3 On, one fifth of the frequency at the resonance 
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peak. The 15th harmonic was included. The transient response 
thus derived is plotted in Fig. 16. It appears to settle down to 
a steady-state value of approximately e. But the actual 
steady-state value is one. The method has failed, because the 
frequency spread is too large. We see from Fig. 15 that the 


frequency spread, as defined in Appendix II, is approximately: 


2 = 

Many other compensating networks were tried in an effort to 
obtain a satisfactory compromise between the three stated 
requirements. In fact, none was found, and it appears to the 
writer that the characteristics of the uncompensated open-loop 
system under consideration make the existence of such a compensator 


impossible. 


In order to continue the investigation of the non-linear 
problem, the actuator, vane and aerodynamic transfer functions were 
arbitrarily changed to make a satisfactory compromise possible. 

The actuator and vane transfer functions were set equal to one, and 


the aerodynamic damping ratio SS was cnanged from .0085 to .%. 


1 
The block diagram of the simplified system is shown in Fig. 17. 
the compensator decided upon was a lead network-amplifier with 


transfer function: 


ae 
=) 

C(s) a lL + Um, 
1+ ie S 

ai 


A Ke, was set at 10 db, again giving a fairly hign value of 


closed-loop peak gain, and requiring that K, = 2.38 db for unity 


B, 


closed-loop steady-state gain. The closed-loop frequency resvonse 


is plotted in Fig. 18. It has a satisfactory frequency spread. 
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SECTION 6: THE VRANSIENT RESPONSE OF THe NON-LINHaAk CLOSED- 
LOOP SYSTEM 
In the simplified system, that is for 
W(s) S V(s) = 
Eqs. 5-1 and 5-7 become, respectively, 


E(t) 


K, X(t) - <(t) (6-1) 


and 


y(t) 


i, CD) att) (6-2) 
Equations 6-1 and 6-2 are applicable to the non-linear system as 
well as to the linear one. 


Using the same notation as in Section 2, the differential 


equation for the aerodynamics is: 


7(D) x = u(D)» - AX? (6-3) 
Where 
(Deere een ts (6-4) 
1 u 
U(D) = K, we (1+ 7D) (6-5) 
dl 


Substituting (6-1) in (6-2) gives 
h(t) = Kc) [ K, X(t) - x(t)] (66) 


Substituting (6-6) in (6-3), 


7(D) X = K, C(D) U(D) | k, x, - «| yee (6-7) 
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This may be rearranged in the form: 
[ 2) + K, c(D) u(p)| x = kK, K, C(D) UD) *_-£,A Se (6-8) 
2 Zao D 


Now an approximation to the solution of kg. 6-8 is the solution 
xX, (t) of: 


| 2(D) + K, C(D) u(D)| ~, = kk, c(D) UD), (6-9) 


In transfer function form, Eq. 6-9 is 


=< ™ K, K, C(s) U(s) 

Xp = Ws) + K, C(s) U(s) 
U(s) 
TS 


Ky K, C(s) 
- 1+K, C(s) us 2 sae 
2 T(s 


Comparing Eqs. 6-4 and 6-5 with Eq. 5-6 we see that 


Us 





Ms) = Ky P(s) (6-1) 
Therefore: 
x \ ae C(s) P(s) 
Os) = 1 ao (6-12 ) 
apart 1+K, K, C(s) P(s) 


The right hand side of (6-12) is simply the closed-loop transfer 
function of the linear system (cf. Eq. 5-9), so that X(t), the 
first approximation to the transient response of the non-linear 

system, is the response of the linear system - not an unexpected 


result. 
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Now by analogy with the argument presented in Section 2, we take 


the second approximation to be: 


X(t) = R(t) (6-13) 
where 
R (%) (6-14) 
(sy 
Ae, is the steady~state output of the non- 


linear system, 


(<X) 


Sq, is the steady-state output of the linear 


system (equal to the magnitude of the yout 


in the present example). 


The third approximation to c(t) is X(t), the solution 
Or $ 


ira ey 3! = B, 
[ 2(0) ao C(D) u(p)|~., = K, kK, C(D) U(D) < (a ~ (6-15 ) 
Or, if St) is the solution of: 
: 3 | 
[ 7(D) + K, C(D) u(p)| Sore a - (3 ma (6-16) 
then by superposition of the solutions of (6-9) and (6-16), 


HX, = XK + &y (6-17) 


The solution of Eq. 6-16 follows the same lines as that of 


Eq. 2-16: we have oo) in the form of a Fourier series: 


A 2 
caSE Go) = Ix,| | * ie s pee: sin (nwt 9 B) | 


where lx, | is the magnitude of the step input, and the A's 
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and B's represent the frequency response of the closed-loop system. 


As before we let 


ers pee 2 FS arate 
= 3 +4 mn opt (nwt + 6B 
so that 
(t) 
sl = |x),| a 


From Eq. 6-13: 


x(t) = R&X)(t) = Rie<,,| x) (t) 
and 
x, = R° Ix]? ae (6-18) 


We find the Fourier series for — : by the method of Section 2, 
and then, by Eq. 6-18, the series ee x, 2 By fq. 6-16, &X 


3x 
is found by applying the frequency Speers 
H(iw) = (2 (6-19) 
Tiw) + K, C(iw) U(iw) 
to the series for ~%,, H(iw)) may be rearranged, using 
Eq. 6-11, to simplify its calculation. It becomes: 
B K, Ky Kz C(iwd) P(i Ww) 





H(iw) = K, K, C(iw) Uliw) 1+kK, K, Cliw) P(iw) (6-20) 


the bracketed term is the linear system closed-loop frequency 


response. 


The next approximation to «(t), if it is desired to go 
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further, is <(t); the solution of: 

Ji j — _ x ? oe 
| °@) + K, c(D) vp) SF K, ke c(D) U(D)*, re ; (Gael) 
Adding and subtracting px, on the right hand side of Eq. 6-21; 


[ 2(0) +k) C(D) u(D)| “~, = K, k, OL Gay IE EDY, aa 


5 3 3 
apo, Bk Hex, (6-22) 
Let S(t) be the solution of: 
3 3 | 
[ 2) + K, C(D) u()f x = on ass = x, ) (6-23) 
Then by superposition of the solutions of (6-15) and (6-23), 


xX, = Xx. + Xs (6-24) 
To solve Eq. 6-23, the Fourier series for ( x.” - x) is 

found by numerical harmonic analysis in exactly the same way as 
outlined in Section 2 and demonstrated in Section 3. ‘the 

frequency response H(iw) is applied to it to give et and 

7 is found by adding the evaluations of x and Xx 


Further iterations follow the same lines. 


Having calculated (+) and %,4(t) for a step function 
of any magnitude, it is easy to find their values for any other 
input magnitude. Let the subscript "a" denote the case which 

has been calculated and "b" the case which is to be calculated. 


The relationship between the linear solutionsis, of course, 


Xx, (4) = Ip, | ee (t) (6-25 ) 
b |X. | a 
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XK, (t) 


R, oX, (t) (6-26) 


and 


NO 
o~ 
oe 
w 
li 


ox = 
: R, 1, {#) (6-27 ) 


o< 
xX, (t) = 2, | Dol ox, (t) (6-28) 


Eliminating &, (t) between (6-26) and (6-28); 


a 
R, | eX 
x, (t) = b | Ds X, (t) (6-29) 
b R, lx, | a 
a 
Thus 
R, [xX 3 
acy 7 » | Do | O53 (6-30) 
b a 
Bee Xp | 


Due to the linearity of Eq. 6-16, 


Axx 


’ X 3% (6-31) 


Re |X | 
a 
Having found XX, (t) and X., (t), X%., (t) follows from 
b "py o 
Kq. 6-17. No new Fourier evaluations have been performed. 


The only simplification possible in the next iteration is 
that the evaluation of x,” follows from Eq. 6-13. x5, 
must be found by et aaa point by point, and a new harmonic 


analysis must be performed. 





47 


We now apply this procedure to the simplified control 


systen, first for a step input %. of magnitude .3 radian. 


D 
The fundamental frequency W., is taken as .38 Hn,» one fifth 


of the closed-loop Be et cae peak in Fig. 18. The 
tabulations of Fourier series and frequency responses for this 
example will be found in Table XIII, and the evaluations in 
Table IXV. The Fourier series for a square wave, magnitude .3, 
is first tabulated, then the closed-loop frequency response. 

The Fourier series for the output of the linear systen, of (t), 
is found by applying the closed-loop frequency response to the 
Square wave. We find the steady-state output of the non-linear 


system from Eq. 5-ll, which, upon substitution of the constants, 


becomes: 

(Mon + 56 (XZ. = [Xp (6-32) 
For 

(5 = |xX,|= 2 
find 

(en = nee! 


Thus, from Eq. 6-14, 


seo 


a ar 


— eee Pei 


x A(t), found by Eq. 6-13, is evaluated. 


The Fourier series for x? , found by the method developed 
v 


in Section 2, is tabulated. 


Rr [& |? = (.957)? (.3)? = .02365 


By Eq. 6-18, ~%? is found from ae and tabulated. The 

frequency response H(iwW) is tabulated, and the fourier series 

for mere” is obtained by Hq. 6-16, tabulated, and evaluated. 

By Ko. 6-17, the evaluation for H(t) is found by adding those 
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for oF and ere 


x,t) and ox, (t) are plotted in Fig. 19. We note that 
x(t) is not very different fron %,(t) and is therefore 
probably a close approximation to o<(t). The effect of the 


non-linearity is small, causing a very slight increase in freouency. 


To illustrate a case where the effect of the non-linearity 
is large, we take (x, | = l radian. Although this is an 
unrealistic value of angle-of-attack, the result will be equivalent 
to that of a missile at low angle-of-attack but with a larger 
non-linearity. It is simply convenient to calculate it in this 
manner since the equations have been formlated on this basis. 
The evaluations are given in Table IXV. From Eq. 6-32 find 


(<<) = .757 


oN 


Using the subscript "b" to denote this case, we have, from 
Eq. 6-14, 


R= stot Biss, 


The subscript "a"! denotes the .3 input case, already evaluated. 
R, <1 

——--— = Avot) Ay = 2.63 

Re |%<p,| (.957) (.3) 

The evaluation of of, ft), found according to Eq. 6-29, is 

tabulated. 


R, leo l\? 
oe I< | 


Xa, 6b) is found by Eq. 6-31, and X 


BON) We we ele 


- O-l/. d 
a by Eq. 6-17 en an 


3, are plotted in Fig. 20. 
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One additional graph is presented in Fig. 21. Here the third 
apvroximations corresponding to both the .5 and the 1.0 input have 
been plotted, as has that for the .65 case, the calculations for 
which are not shown. These solutions have been "normalised" by 
dividing each ordinate by the corresponding steady-state value 
(X) %,(£)/() , which, of course, is the same for all 
inputs, has also been plotted. This is a useful method for 
visualising the effect of increasing input magnitude on ~<,(t). 

It would also be useful for plotting other approximations, or the 


exact solutions. 


Improved approximations will not be demonstrated here, and 
the determination of the number of approximations necessary for 
acceptable convergence remains for the worker who thinks this method 


might be useful. 


Finally, it 1s interesting to note that a linear closed-loop 
system with two inputs can be described which is equivalent to the 
non-linear system. The block diagram of this system is shown in 
Fig. 22. Comparison with Eq. 6-8 will prove its equivalence. 

Fig. 22 is, in fact, the general outline of a possible analog 
computer set-up where ax? is found by putting c& through a cubing 


device. 
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DISCUSSION 


It has been mentioned that a quick method for determining the 
effect of non-linear stiffness was desired. specifically, it was 
hoped that few terms would be required in the first iteration, and 
that one iteration would be sufficient for acceptable accuracy, at 
least in the first overshoot region. In fact, the number of 
Significant terms was not small, and a single iteration was not 
sufficient no matter how many terms were retained, except for small 


inputs. 


The frequency spread requirement, discussed in Appendix II, 
places a limitation on the range of systems to which the Wass, 
Hayman method may be applied. And the Wass, Hayman method must be 
applicable for the present method to work. If the frequency spread 
is of the order of 20, the number of terms in the Fourier series 
could be increased to handle it, but if the original control system 
designed here is at all typical, the frequency spread is likely to 
be at least 50, requiring at least 25 terms in the Fourier series 
for the linear response, and more for the non-linear response. 

The need would have to be great to justify the effort of evaluating 


such series. 


In Section 5 an unrealistic value of aerodynamic damping was 
arbitrarily chosen to produce a control system to which the present 
method could be applied. It is possible that a more justifiable 
change in tne control system could nave been made, but this was not 
done due to a necessity to find an acceptable system guickly. It 
remains to be seen whether the present method can be used in any 


practical cases. 


Although only "hard" stiffness has been considered, i.e. 
positive (3 » the method is equally applicable to soft stiffness, 
or negative (3 : In the latter case the presence of static 


instability at high angle-of-attack must be considered. 
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CONCLUSIONS 


l. Use of the iterative technique developed is 
capable of giving the transient response, to an acceptable degree 
of accuracy, if the frequency response of the system meets the 


frequency spread requirement. 


EL In an example, acceptable accuracy was obtained 
after two iterations. The number of iterations necessary when 


the non-linear effect is large was not determined. 


On The labor required to evaluate the Fourier series 
and to perform the numerical harmonic analyses necessary in the 
second and subsequent iterations is considerable, thus limiting 


the practical value of the method. 
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APPENDIX I 


CALCULATION OF THE ABRODYNAMIC CONSTANTS 


Brown-Edwards, in Ref. 1, found, for the missile under 
consideration and for a c.g. position 8 body diameters aft of the 


body apex, the following values: 
C = ZV 2 (leap) 


c ane ee eG (1-2) 
are | 


Fixing the body diameter 1 at one foot, he estimated for a 


typical missile: 


2 
lb. sec 
m = TA 05 rt. 
and 
is "2450 eee 
2 


The reference area S if QO ft. 


We take the case of flight at 50,000' at ach 3. she air 
density ( is .000362 lb. sec*/tt* and the speed of sound c 
is 968 ft/sec. Thus the velocity is: 


V= 3 (968) 2900 ft/sec. 


The dynamic pressure is 


4(.000362)(2900)* = 1520 = 


&> 
I 
n]— 
re 
bh 
it 


and 


QS = (1520)(9) = 13,700 lb. 








oc, 
Q S eee (13,700)(2.87) = — 39,406 1b. 


ca¥) 


oe 
Ox 


Z.. = 


where OC /ax is obtained from Eq. I-l. 


oho | 
zZ, Ate oy = hor + = = (13,700)(- .222) = 3,040 1b. 
where OC./On is obtained from Eq. I-l 
MM. = OM = as 
h an 7 Qo = (13,700)(1)(.97) = 13,300 1b.ft. 


where OC, /On is obtained from kq. I-e. ‘he moment due 
Cee: 


to « is, from Eq. I-2, 


1" (oX) Oat eine ace x) 


(13,700)(1)(- 114K - 2.650) 


(-15,600 & - 36,306%°) 1b.ft. 


Comparing this with Eq. 1-14, we see that 


M x =i OOou baat. 


and 


M 


5 = 36,300 lb.ft. 
XX 


Using Eqs. 1-17 to l-edl, 


—2 
Worse =e = Mex = 15,600 = 67.8 sec . 


I 250 











oe — top eae) ge 
3, = -_&« = __ 39,400 
a" me 2(12.05)(2900) (8.23) 
= .0685 
K, = My 2 13,300 Es: 
lrg ee 
” (230) (67.8) 
. =A ee 
C aad 2 a 
mV K, Un (12.05 )(2900)(.853) (67.8) 
1 
= 001505 sec. 
oe 36 300 -2 
(3 = - —— = Sa 158 sec 














APPENDIX II 


THE WASS, HaYMAN METHOD OF DERIVING THE TRaNSIBNT RESPONSE OF 


A LINEAR SYSTEM FROM THE FREQUENCY RESPONSE 


Reference 2 presents a method of deriving the approximate 
response of a linear system to a step input from the frequency 
response. The method is to find the transient response to a 
Square wave input instead of to a step input. The justification 
for such a manouver is best explained with reference to Fig. 23, 
Figure 23a is the graph of a step input, and Fig. 25b is the 
corresponding output of a typical system. Figure 23c is a square 
wave input of period T, and Fig. 25d is the corresponding output. 
It is evident that if the half-period T/2 is greater than the 
time required for the transient to completely damp out, then any 
of the half cycles of Fig. 23d is the same as the transient part 
Grebe. 250, Of course in linear systems the transient never 
completely damps out. Therefore Fig. 23d is really an 
approximation to Fig. 23b, its accuracy depending on the extent 


to which the transient subsides in the time 1/2. 


The reason for using a Square wave is that it can be expressed 
as a fourier series. The Fourier series for a square wave F(t) 
with equal mark-to-space ratio, minimum value zero, maximum value 


unity, and fundamental frequency We » can be written: 
1 2 
F.(t) = oft _ > 1 sin nwt 
al a n F 
aloe 25 on eoeey a] 


(Wass and Hayman use the symbol n differently - this will be 
explained shortly). If F(t) is the input to a linear system 
the magnitude and phase of each component will be changed 


according to the frequency response of the system, and the 








output will be 


A A 
Huey) O > a 
O = sas eo oie (nw a B) 


A Iro 


n=. or oe eee C8 


where AL is the gain of the system at frequency nw, 


BO is the phase shift of the system at frequency nWd,.. 


In order to use this method numerically, the summation mst 
be limited to a reasonably low number of terms. Wass and Hayman 
state that it has been found by experience that a satisfactory 
compromise between labor and accuracy 1S acnieved by considering 
up to the llth, or sometimes 15th, harmonic of the fundamental 
square wave frequency. Having decided on the number of terms 
to be summed, a value for the fundamental frequency Ws must be 
chosen. Wass and dayman suygest that for systems whose gain 
curves exhibit a definite peak, WW. should be one fifth of the 
frequency of the peak. suggestions are also made with respect 


to other types of systems. 


wass and Hayman do not point out that taking the summations 
up to even the 15th harmonic will fail to give satisfactory 
accuracy in certain types of systems. They do infer that the 
transient response is mainly determined by the characteristics of 
the frequency response at frequencies lower than that at which 
the gain is 15 db below the zero frequency gain. They do not 
point out that the gain and phase shift at the fundamental 
frequency must not differ greatly from their zero frequency values. 
for conciseness we may use the term "frequency spread" to denote 
the ratio of the -15 db frequency to the frequency at which the 
gain or phase first differs appreciably from its zero frequency 


value. Thus the Wass, Hayman method, summing up to the 15th 








harmonic, is not applicable to systems whose frequency spread is 


greater than 13. 


An important part of the Wass, Hayman report is a table, the 
"we" Table, which may be used to simplify the evaluation of Fourier 
series. This is a table giving the values of nwt! in degrees, 
for 59 values of non-dimensional time t', and for odd values of 
n from 1 tol13. Even values of n are not included because the 
Fourier series for the responses of linear systems to square wave 
inputs contain only odd harmonics. Since the evaluation of even 
harmonics is required in this thesis, an extension to the Y Table 
is required, and this is given, for even values of n from 2 to 8, 
in Table XV. Wass and Hayman's use of the symbol n precludes 


the possibility of even harmonics, and this is the reason it is 


used differently here. 


Wass and Hayman present a design for a circular computer which 
reduces considerably the time required to evaluate a Fourier 
series. The computer suggested consists of a circular base with 
various lines and scales marked on it, two arms which pivot at the 
center of the circular base, and a cursor which slides along one 
of the arms. The author constructed and used such a computer, and 
eventually devised a different design which was superior in several 
respects. Fig. 24 is a photograph of this computer. It was 
originally a "Wind and Navigation Calculator Mk. I", and all the 
scales around the periphery are superfluous except tiie 360 degree 
one. It consists of an opaque base and a transparent top, 
pivoted at the center. The horizontal lines, the numbers running 
down the center, and the arrowhead are etched on the base. The 
360° scale is etched on the top. The top is roughened so that it 


will take pencil marks. To evaluate the term 
A = B sin (nwt! + ©) 


the procedure is as follows: 








1. Rotate the top until the arrowhead points to 
(90° = §@) on the 360° scale. 


2. Make a pencil mark, on the top, at the point 
where the value of B appears on one of the vertical scales 


above the center. 


3. Rotate the top until the arrowhead points to 


nwt! on the 360° scale. 


a Follow the new position of the pencil mark along 
the horizontal lines to the same scale used to place the mark, 


and read A. It is positive if above center, negative if below. 


5. Repeat steps 3 and 4 for other values of nwt. 


The new computer 1s superior to the Wass, Hayman version in 


these respects: 


a. The two arms of the Wass, Heyman computer must be turned 
through several revolutions without moving relative to each other, 
yet tne angle between them must be adjustable. A design 
incorporating such a feature may be difficult to achieve, particularly 


for the home builder. The new computer obviates this difficulty. 


b. The angle& between the arms and the position of the 
cursor must be reset for each new term of the Fourier series. 
With the new computer a mark for each term of the Series can be 
made at the outset, each one being labeled appropriately. ‘This 
may not seem important, but in practice it is. Mistakes can be 
corrected easily because the setting is not "erased" for the next 
term as with the Wass, Hayman version, and evaluations can be made 
for additional values of time, without resetting, after several 


points have been plotted. 








APPENDIX III 


GRAPHICAL ADDITION OF LIKE HARMONICS 
If, in the expression 
e sin (wt + ©) +f sin (wt + ©) +g sin (wt +y) 


e, f, g, ©, P, and yw are known constants, then an equivalent 
expression is 


Wh ciao n ) 


where h and nh are determined according to the construction 


below. 





This method is valid for any number of terms in the original 


expression. 





APPENDIX IV 


TRANSFORMATION OF THE CUBIC EQUATION 


Given the equation: 


2 d 
oes + ay - + aX + ee eee (I~1) 


Make the substitutions: 








1 
oe a ¢ 
O 
and 
6 
ni Nou *( 
Then Eq. I-l becomes: 
dy ~ 4 dy ty ty" 2 ae 
Eee J dc 6 
We shall call the quantity 
9 - ate 
O = maser.) 4 
ate 
O 
the non-dimensional input. The coefficient = is equal to 
O 


twice the damping ratio §. 











APPENDIX V 


DESIGN OF THe INCIDENCE VANE 


The vane is conical, and it is necessary that the c.g. be as 
far forward as possible. The configuration to be used is that 
shown in Fig. 25. The nose will be solid steel, and the base will 
be a thin shell of light-weight plastic. The optimum ratio of 
nose length to tail length, for most forward c.g. position, is that 
which places the c.g. at the junction of the nose and tail. The 


equation for this condition was found to be: 


. i a sin © 
say (2 oe 3 oe n) 


tifct 


where is the base thickness, 
is the cone length, 


is the density of the nose material, 


is the cone semi-angle, 


1 

L 

(i 

Po is the density of the base material, 

iS 

= is the distance of the c.g. from the apex. 
The density of steel is about .3 lb/c in., and that of a high- 
temperature silicone plastic is about .05 1b/cu in. (p. 66 of 


Saar. Using these values, the compromise chosen was: 


Sees oe 
n= 3 
oe. 
t= 0614" 


It is questionable whether the base is thick enough for adequate 


strength, but it will be assumed that such a design is practicable. 








APPENDIX VI 


CALCULATION OF THe VANE CONSTANTS 


oe 2 a 


Ref. 5 gives experimental pitching moment data for cones in 
supersonic flow in the form of graphs of C,, versus XX, where 


the pitching moment about the nose is 


a 
Uvesne = SS, 1 Ge 
and 
ie eS base radius, 
l = cone length. 


The curves of Cr versus & were virtually linear for the conditions 
applicable to tne present case, and the value of 9¢,/ Ox about 
the nose was .020. the center of pressure remained very close to a 
position 2/3 1. back from the nose, so that the moment curve slope 


for a cone hinged at 1/3 1, would be half of .020, or .Q10. 


From A,pendix I: 


q 1520 1b/ft* 


rom Apoendix Vv 


L. =. “25: tts 
Re = J0107 £20. 
thus: (2%) ; (2x | 
G = - ee = ~ 11TQ2 aE ==) 
Ht aK b Cc xX . 
2/3 e 1/4 


~ 77 (1520)(.0107)(.25)(.01) 


ll 


ae nize Tet aeleoe 


1.70 x 107° ft.lb.sec* 


L. 


(“ound by integration, calculation not shown). 


uw) = - a = ieee = 7,530 seo“ 
se ~5 
ly 1.70 x 10 
| -1 
WwW = 86.9 sec 
n 
2 
5, = 5 


(Chosen arbitrarily on the assumption that the vane damping device 


can be so designed. ) 
C, = 2 ty M252 = 2(1.70 x 10-)(86.9)(.5) 


= .00148 fi.1b.sec. 
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